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We consider a hole-doped semiconductor with a sharp boundary and study the boundary spin 
accumulation in response to a charge current. First, we solve exactly a single-particle quantum 
mechanics problem described by the isotropic Luttinger model in half-space and construct an or- 
thonormal basis for the corresponding Hamiltonian. It is shown that the complete basis includes 
two types of eigenstates. The first class of states contains conventional incident and reflected waves, 
while the other class includes localized surface states. Second, we consider a many-body system in 
the presence of a charge current flowing parallel to the boundary. It is shown that the localized 
states contribute to spin accumulation near the surface. We also show that the spin density exhibits 
current-induced Friedel oscillations with three different periods determined by the Fermi momenta 
of the light and heavy holes. We find an exact asymptotic expression for the Friedel oscillations far 
from the boundary. We also calculate numerically the spin density profile and compute the total 
spin accumulation, which is defined as the integral of the spin density in the direction perpendic- 
ular to the boundary. The total spin accumulation is shown to fit very well the simple formula 
5tot <x (1 — m_t/m_r/) 2 , where m_L and ran are the light- and heavy-hole masses. The effects of 
disorder are discussed. We estimate the spin relaxation time in the Luttinger model and argue that 
spin physics cannot be described within the diffusion approximation. 

PACS numbers: 

I. INTRODUCTION 

Hole-doped semiconductors are a very well studied and industry developed class of materials. The fundamental 
description of these materials is usually based on effective models such as the Kane model or Luttinger models 
which capture most of the properties of a semiconductor. A key ingredient of these models is spin-orbit interaction, 
which couples the momentum with the orbital and spin degrees of freedom. It should be noted that the latter 
degrees of freedom in semiconductor systems attracted attention only very recently, when it was recognized that 
the spin-orbit coupling may lead to the possibility of spin control by electric means On one hand, the predicted 
spin-charge coupling opens a possibility for new useful spintronics applications^^ On the other hand, it leads to a 
variety of new theoretical problems, which need to be clarified in order to understand the relevant experiments^^. 
The theoretical description of the intrinsic spin-Hall effect j&sLiSi which is one of the spin-charge coupling phenomena, 
relies on an elegant mathematical structure known as the Fermi surface Berry's phase .iLi^iiii This structure originates 
from the spin-orbit splitting of the bands. Band crossings become sources of a fictitious magnetic field in momentum 
space, which leads to a non-trivial Berry's phase and may affect certain observables. In particular it leads to an 
anomalous contribution to the velocity operator. Consequently, any quantum mechanical operator, which involves 
the velocity acquires an anomalous contribution. An important example is the spin current operator (usually defined 
as a symmetrized product of the spin density and the velocity), which may acquire an anomalous component as 
well, if an electric current is present lii^i^iiSiiLiSiiS The existence of the anomalous contribution to the spin current 
perpendicular to the electric current is an important prediction of the spin Hall effect theory. However, a direct 
experimental check of this prediction, while possible in principle, is not straightforward because the spin is not 
conserved, and as such, the spin current has no obvious physical meaning^SiSI Alternatively one can experimentally 
probe current-induced equilibrium spin density and that is what usually is measured in experiment. It is therefore 
desirable to develop a theory, which would allow one to calculate observables such as boundary spin accumulation 
directly and provide a clear understanding of the physical processes behind this effect. It is also desirable to search 
for other possible manifestations of the Fermi surface Berry's phase apart from the anomalous velocity. 

In this paper we consider a three-dimensional hole-doped semiconductor described by the isotropic Luttinger model 
and in the presence of a boundary. We mostly discuss the clean limit when no impurities are present or, alternatively, a 
disordered system but only at small (ballistic) length-scales. We note here that the application of an electric field to a 
perfectly clean system would result in a non-equilibrium state and time-dependent spin density. To access equilibrium 
spin-Hall physics we assume that either the voltage drop occurs only in the contacts, or that there is a relaxation 
provided by impurities. In both cases there is an equilibrium charge current, which determines spin accumulation 
near the boundary. 

Evaluation of the current-induced spin density involves solving the single-particle Schrodinger equation for the 
Luttinger Hamiltonian in a half-space. The corresponding solution contains a few important features, which are quite 
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different from the usual single-band quantum-mechanical problem. While the states of the bulk Luttingcr model can 
be classified by quantum numbers corresponding to the double-degenerate heavy and light-hole bands, the boundary 
does not conserve these quantum numbers and mixes up different bands. For example, if a heavy-hole with positive 
chirality propagates towards the boundary in the direction close to normal, it gets reflected in all bands, which are 
heavy- and light-holes with positive and negative chiralities, and the reflected light holes have the angle of reflection 
(measured from the normal to the surface) larger than the angle of incidence. An important property of the solution 
is that for large enough angles of incidence, light holes arc not reflected from the boundary at all, but instead get 
localized near the surface. These states are similar to the Tamm states^ which appear in crystals due to an abrupt 
change of the electronic band structure at the boundary. We also point out that one can draw an intuitive analogy 
with optics by imagining that the heavy holes occupy a half-space with high index of refraction, while the light holes 
occupy the other half with a lower index of refraction. A wave propagating say from the medium with the high 
index of refraction (heavy hole) may get reflected from the interface (remains a heavy hole) or may get refracted 
and propagate in the other medium (becomes a light hole). For large enough angles of incidence, one expects total 
internal reflection, which is somewhat similar to the appearance of localized light-hole states in our language. It can 
be shown that these localized light-hole states contribute to spin accumulation if a current is present (below we study 
accumulation of the total orbital momentum, but occasionally call it "spin accumulation" for the sake of brevity). 

To qualitatively understand the physics of the boundary spin accumulation in a many-particle system it is useful 
to recall the well-known problem of a free Fermi gas in the presence of a boundary. The boundary (or any other 
perturbation for this matter) leads to the Friedel oscillations in the particle density with the period of two Fermi 
momenta. We note that the integral of the density over distance reduces to the bulk density, i.e., there is no boundary 
charge accumulation since the latter is conserved. A similar oscillatory behavior of the spin density may occur in 
spin-orbit coupled systems (e.g., described by the Luttinger model) if a current is present. There arc, however, two 
important differences: (i) The existence of multiple periods of Friedel oscillations (or beatings) due to two distinct 
Fermi momenta corresponding to the light and heavy holes; (ii) The non-zero integral of the spin density due to the 
non-conservation of spin. Summarizing the above qualitative discussion, we argue that non-zero spin density, which 
appears near the boundary in response to an applied current can be viewed as current-induced Friedel oscillations 
with non-zero spin accumulation originating from the localized (Tamm-like) surface states. 

Our paper is structured as follows: In Sec. [H] we solve the quantum-mechanical problem of a particle described 
by the Luttinger Hamiltonian in a half-space. We show that for the angles of incidence greater than a critical angle 
8 C , the eigenstates contain modes that are localized in the direction normal to the boundary. We construct a full 
orthonormal basis for the problem. Even though the formulation of the problem is very straightforward, its solution 
is technically challenging due to a cumbersome matrix structure of the Hamiltonian. Subsections III Al and III Bl arc 
devoted to a formal proof that certain symmetrized combinations of incident and reflected waves constitute a full 
orthonormal basis for the Hamiltonian. A reader not interested in the formal proof should skip to Sec. II I II 

In Sec. IIIII we address the question of boundary spin accumulation in response to an external current. Using the 
orthonormal basis constructed in Sec. El we obtain numerically the spin density profile near the boundary for various 
values of the ratio of the light and heavy hole masses, £ = m^/TO^. We also extract analytically the large-distance 
asymptotic behavior of the spin density, which is shown to oscillate with three distinct periods 2k^\ 2k^\ and 

+ k^ (where k^ and kp are Fermi momenta of the heavy and light holes correspondingly) and decay away 
from the boundary as a power law, cx 1/r 2 . In subsection llll Bl we define a quantity, which we call spin accumulation, 
by integrating the spin density over distance in the direction perpendicular to the boundary. We obtain the spin 
accumulation numerically and show that the total spin accumulation fits the simple formula Stot oc (1 — £) 2 . 

In Sec. IIVI we discuss qualitatively the effects of disorder on the current-induced Friedel oscillations. Using the 
analogy between the spin accumulation effect and the usual Friedel (or RKKY) oscillations in disordered metals, we 
argue that one should expect an interesting behavior of the current-induced Friedel oscillations in the spin density 
in a disordered system. We argue that while the system-wide average value of the boundary spin density decays 
exponentially away from the boundary as e~ r l l (where / is the mean free path), the higher moments and the typical 
spin density still decays as a power law r -2 . However, the latter has a random sign and the spin accumulation (i.e., 
spin density averaged over large enough distances) decays exponentially as e~ r l Ls (where L s is a spin relaxation 
length). We calculate the spin relaxation time and show that for all reasonable values of the spin-orbit coupling, the 
spin relaxation time is very short and has a universal value t s = 3r/2. This short spin relaxation time implies that the 
spin relaxation length is of the order of the mean free path and, therefore, the hydrodynamic diffusion approximation 
docs not apply for the Luttinger model. 
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II. EIGENSTATES FOR THE LUTTINGER HAMILTONIAN IN HALF-SPACE 



The physics of a hole-doped semiconductor with diamond or zinc-blende structure is often described by the effective 
Luttingcr Hamiltoniani: 



T^Lut — 7; — 
2m 



(i) 



where m is the effective mass, k = (k x ,k y ,k z ) is the momentum and S = (S x , S y , S z ) represents the total angular 
momentum 3/2 of the atomic orbital, i.e. the sum of the orbital angular momentum and the spin. The total angular 
momentum can be represented by three 4x4 matrices with explicit expressions given in Appendix [S] For simplicity, 
we consider the spherically symmetric model described by one Luttinger parameter^, 7, and we choose the units so 
that fi = 1. 

For a system with translational symmetry, the Hamiltonian Q can be diagonalized in a basis in which the helicity 
operator A = k • S/fc is also diagonal. For a given wave vector k, 7init has two double degenerate eigenvalues 

ejf ( k ) = L^l k ^ = JL for A = ±-, 

2m 2mn 2 

e L (k) = l -±^lk^^- for A = 4- (2) 

2m 2m l 2 

These two bands are referred to as heavy-holes and light-holes, respectively The corresponding eigenfunctions can 
be expressed in terms of a four-component spinor, 

$ kA (r) = ^=e ikr £AK), (3) 



where f2 is the volume of the system, the label A is (±3/2) = (H±) for the heavy-holes and (±1/2) = (I/±) for the 
light-holes, and the spinors U\(n^) depend on the orientation of the wave vector, n k = k/|k|. Their explicit forms 
are given in Appendix lAl 

Next we consider a similar problem for a system defined in half-space, i.e. in the presence of a sharp boundary in 
the z direction defined by the potential 

nr) = {° * Z >° (4) 
w I 00 if z < ' 

Solving the quantum problem described by the Luttingcr Hamiltonian Q in the presence of the boundary is still 
a straightforward exercise. However, finding an orthonormal basis for this system, represents a technical challenge 
because of two main reasons: 1) The new eigenstates are, in general, linear combinations of bulk eigenfunctions < E > kA( r ) 
and are not automatically orthogonal on each other, and 2) The proper counting of modes is difficult, as the standard 
technique of imposing periodic boundary conditions and discretizing the momentum is not immediately applicable for 
eigenstates involving combinations of heavy and light holes. In this Section, we will direct our effort toward solving 
these difficulties with the goal of constructing an orthonormal basis for the Luttingcr Hamiltonian in a half space. 
Let us notice that while k x and k y are still good quantum numbers for the full Hamiltonian TL = H-Lut + ^( r )j k z and 
the helicity A arc not. The wave functions have to vanish at the boundary, i.e. for z — 0, while for positive values of 
z they are superpositions of incident and reflected waves. Let us consider the case of a A ~ +3/2 incident heavy-hole 
(see Fig. QJ. We can always choose the coordinate system so that k y = and k x > by properly rotating the axes. 
The full wave function corresponding to the reflection process represented schematically in Fig. ^ has the following 
form for positive values of z, 

lbi H+) (r) = t^L. {U H+ (0y k * z + A 1 U H+ (n - 0)e- ik ' z + A 2 U H -{* - 0)e^ z 



+ B 1 U L+ (n - 4>)e^ z + B 2 U L -{ir - 0)e^ 2 } , (5) 

where C is a normalization factor. All the reflection coefficients Aj(0) and Bi{9) are in general non-zero and are 
determined by the boundary condition ipk(z = 0) = 0. The parameters that describe the reflected waves can be 
determined from the parameters of the incident wave using the momentum and energy conservation laws. Explicitly 
we have for the heavy-holes 

V - k ■ 

k' y = k y = 0; -> e' = 7r-e, (6) 
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FIG. 1: (Color online) Scattering of an incident heavy-hole plane wave. By convention, the coordinate system is chosen so that 
k y — and k x > 0. The angles 8, 6' and <p that describe the direction of propagation for a given wave are defined as the angles 
between the z-axis and — p, where p 6 {k, k', q'}, and belong to the interval [0, n]. The regime represented here corresponds 
to 9 < 9 C (see text) and the incident and reflected waves are described by the wave-functions defined by Eq. ©. 




FIG. 2: (Color online) Scattering of an incident heavy-hole plane wave with the incident angle 9 > 9 C - The reflected light-holes 
propagate parallel to the boundary and are localized in the z-direction. The wave function describing the reflection process is 
given by Eq. Notice that the localized modes are described by two independent spinors Vi(x) and Vsj(x) ( see Appendix IXt . 



where 9 = arccos(— k z /k) <E [0, 7r/2] is the angle of the incident heavy-hole. Similarly we obtain for the light-holes 




0; -> 0' = 7T - = 7T - arccos I Jl-^M I , 

sm 2 (9)} 1 / 2 ; £ 

where £ = mz,/mij = (1 — 2j)/(l + 27) is the ratio between the light-hole and heavy-hole masses. Consequently, 
the wave function (JSJ) is an eigenfunction of the full Hamiltonian with an eigenvalue = fc 2 /2m# = q 2 jlrtii,. This 
solution exists as long as q z is a real number, i.e. for incident angles smaller that the critical angle 

9 C = arcsin(V^). (8) 

For 9 > 9 C there arc no light-holes propagating in the z-direction but, instead, the scattering problem has solutions 
that are localized at the boundary. This situation is illustrated schematically in Fig. [2J The localized modes are 
characterized by an imaginary wave vector q z = iQ and the corresponding wave function becomes 

^ + \t) = ^{U h+ {6) e^ 2 + A 1 U H+ (n-9) e~ ik ' z + A 2 U H -(n - 9) e~^ z 
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FIG. 3: (Color online) Family of reflection processes characterized by a given pair of reflection angles (9, <f>{9)). The color code 
is the same as in Fig. Panel A represents schematically the same event as Fig. the scattering of a heavy-hole with helicity 
A = +3/2, while panel B represents the scattering of an incident heavy-hole with opposite helicity, A = —3/2. Panels C and D 
illustrate "re-combination" processes where superpositions of heavy-holes and light-holes are reflected into a single heavy-hole 
mode with helicity +3/2 and —3/2, respectively. For incident heavy-hole angles 9 > 9 C , the light-holes propagate parallel to 
the boundary (<f> = tt/2), while they are localized in the z-direction, and are described by the spinors Vi(x) given by Eq. 1X51 . 



B 1 V 1 ( X ) e- Qz + B 2 V 2 ( X ) e- Qz } 



(9) 



where Vi(x) are spinors describing the evanescent modes and the reflection coefficients are determined, as before, by 
the boundary condition. Explicit expressions for the V*(x) spinors are given in Appendix lAl The angle x an d the 
wave vector that characterize the localized states arc determined again using the conservations laws and we have 



- ky — 0; 
q' z =iQ = ik[sin 2 {9) -£] 



1/2. 



X = arccos 



sin(0) 



(10) 



Notice that all the modes contained in the superpositions that define the wave-function in Eq. JSJ and Eq. J3J arc 
eigenfunctions of the Luttinger Hamiltonian corresponding to the same eigenvalue, 



HLut [U H ±(0) 



e k [U H ±{6) e 



zkrl 



u L ±{4>) (>*(*.*+«.*) 



H hut \Vi(4>) e 



Vi{<f>) e 



ik x x — Qz) 



where 



k x + k z 
2m H 



kl+ql 
2m L ' 

2mi 



if < 9 C 
if 0>0 r 



(11) 



(12) 



Also notice that all the scattering angles and the wave vectors can be expressed uniquely in terms of and k = 
+ ky + kl) 1 / 2 . In addition, the scattering coefficients A, and Bi in Eq. © and Eq. (JHJ are functions of the angle 
of the incident heavy- hole. Their explicit form is given in Appendix iBl 

The wave function (JSJ is not the only eigenstate of the Hamiltonian H = 7iLut + ^( r ) with the eigenvalue = 
k 2 /2rriH that can be expressed as a superposition of propagating modes characterized by the scattering angles (9, 
4>{ffj). For example, a similar state can be obtained by starting with an incident heavy-hole with helicity A = —3/2, 
(HH-). Moreover, we can imagine "re-combination" events, such as those represented schematically in Fig. [3] (panels 
C and D), or the corresponding scattering processes of light- holes shown in Fig. 0] From this family of states one 



6 



can extract sets of four linearly independent cigenfunctions, while all the others states can be expressed as linear 
combinations of these vectors. Our task is to identify such a complete set of orthogonal wave functions that would 
allow us to construct an orthonormal basis for the full quantum problem. The basic idea is to construct this set by 
taking certain symmetric linear combinations of the cigenfunctions represented in Fig. [21 The general expression for 
such a symmetrized eigenstate is 



1 



*£(r) = -=e ik " {[c? U H+ (8) + <% U H .(6)} + [eg U L+ {4>) + c% U L ^)\ e l 



+ [ct U h+ (tt - 0) + eg U H -{it - 8)} e~ tk * z + [c? U L+ (n - 0) + c£ C/ l _(tt - 0)] e^"} 



(13) 



-i+4l 



i 6 {1,2,3,4}, i.e. the weight of each mode is the same for both the incident and the 
reflected waves. In constructing the symmetrized eigenfunctions, we use the fact that the reflection coefficients for 
the family of states described in Fig.[3]can be expressed in terms of the reflection coefficients of the state %j)^ H+ ^ = tp A 
given by Eq. (|B1|I . We write the corresponding wave- functions in the most general form given by Eq. (|13fl . For 
example, the wave- function ip A , corresponding to the process represented in Fig. [SJ'V, will have the coefficients c A = 1, 
c 2 = c f = c 4 = 0; c A = Ax, etc., up to an overall normalization factor that we omit. In Table[I]we summarize the 
relations between the coefficients of the wave-functions corresponding to the reflection processes represented in Fig. 
[3] omitting an overall normalization factor. 
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TABLE I: Coefficients for the family of wave-functions represented in Fig. [3] The wave functions are expressed in the general 
form given by Eq. 11131 . The coefficients corresponding to the reflections from panels B, C and D (see Fig. El are expressed in 
terms of the coefficients for the process A which are given explicitly by Eq. <BH . 

Next we introduce two pairs of symmetrized states defined as 



1 



l 



v/T 



ab 



b 



AC± 



1> 



AC± 



1 



(14) 



where ip AC± = l/VN(ip A ± if; c ) and ip BD T = l/^/N 1 ^ 3 T^ D ), with N and N' being normalization factors that 
insure (tp a \tfj a ) = 1 for both linear combinations. The parameter 'b' is set equal to the scalar product of the two sets 
of linear combinations, b = (ip AC± \ip BD T) , insuring the orthonormality of the eigenstates, 



(15) 



where Sij is the Kronecker 5-symbol. By varying the parameter 'a' we can change the weight of the heavy-hole 
and light-hole modes that contribute to a particular eigenfunction. For a state given by Eq. l|13f) the weight of the 
heavy-holes is defined as 



W, 



HH 



Icil 2 



(16) 



while for the light-holes we have W^ H = 1 — W§ H . We choose the parameter 'a' to maximize the weight of the heavy- 
holes in the eigenstates , I ,1± and the weight of the light-holes in 5 ,2± . The convenience of this choice will become 
clear once we address the problem of counting the eigenstates. The explicit value of the parameter 'a' satisfying this 
condition is 



1 



2W 



IF, 



AC± 



HH 



(17) 



where W = c ac± c BDt + c ac± c BDt + c ac± c BDt + c ac± c BDt . Notice that, while in general the eigenstates 
represent a superposition of both heavy-hole and light-hole modes, in the particular case of normal incidence, 8 = 0, 



FIG. 4: (Color online) Light-hole reflections that belong to the same family as the processes represented schematically in 
Fig. [3] These reflections involve a single incident (panels A' and B') or reflected (panels C and D') light-hole mode. The 
corresponding wave-functions can be written as linear combinations of wave-functions describing the heavy-hole scattering 
processes represented in Fig. |3] Same color code as in Fig. blue - (HH+), green - (HH-), red - (LH+), yellow - (LH-). 
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TABLE II: Relations between the coefficients of the symmetric eigenstates for an incident angle for the heavy-holes 9 < 6 C . 
The wave functions for the pairs of states ** (2)+ and ** (2) ~ (see Eq. O) are written in the general form given by Eq. ltT5t . 
The existence of two pairs of modes labeled by 1 and 2 is reminiscent of the heavy- and light-holes of the bulk, while the degree 
of freedom labeled by '+' and '-' is the analogous of helicity. Notice that all the coefficients can be expressed in terms of two 
sets of four parameters, one set for each pair of modes a = 1± and a = 2±, respectively. 

^E' 1± contains only heavy-holes and ^ 2± contains only light-holes. In addition, for 9 = 9 C , ^ 2± reduces again to a 
superposition of light-holes propagating parallel to the boundary, i.e., 4>(9 C ) = n/2. There is a symmetry between the 
modes contributing to and the states with opposite helicity contributing to . This symmetry translates into 
a set of relation between the corresponding coefficients that is summarized in Tabic [H] 

So far we have addressed the problem of finding a set of symmetric eigenstates in the case when the incident angle 
for the heavy-holes is smaller than the critical angle, 9 < 9 C . Let us focus now on the case 9 > 9 C , when the light-holes 
are localized near the boundary. There are two main differences between the two cases. The first one concerns the 
number of linearly independent eigenstates. Far from the boundary, the system should be locally equivalent to an 
infinite system and, consequently, the number of degrees of freedom should be the same. The set of four independent 
eigenstates that we obtained for 9 < 9 C corresponds to the double degenerate heavy-holes and light-holes of the infinite 
system. In contrast, for 9 > 9 C only the heavy-holes modes propagate in the z-direction and, therefore, we expect 
to have only two independent cigenfunction for a given incident angle. The second difference is a formal one and 
consist in the spinors Vi(x) being complex, in contrast with U\ which are real. Consequently, the coefficients of the 
eigenstates involving localized modes will be complex numbers. The coefficients Ai, Bi for the wave function 10 are 
given explicitly in appendix^ The symmetric eigenstates can be written again as superposition of reflection processes 
analogous to those in Fig. except that this time the light-holes are localized modes confined near the boundary. 
The most general form of such an eigenstate can be written as 

*£(r) = J=e^* {[cl U H+ (9) + c a 2 U H .{6)] e ik ** + [cf V^x) + c% V 2 { X )\ e~ Qz (18) 



8 



[cf U H+ (ir-6) + c% U H -(tt 



-ik z z 



}• 



where Q(k, 9) and x{8) are given by Eq. (|1U[) . Again, the coefficients corresponding to a family of reflection processes 
similar to that represented in Fig. |3| are not all independent and we summarize the relations between them in Tabic 
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TABLE III: Coefficients for a family of wave-functions with localized light-hole modes analogous to the scattering processes 
represented in Fig. The wave functions are expressed in the general form given by Eq. I|18|l . The complex coefficients 
corresponding to the reflection of a heavy-hole with helicity -3/2 (analog to the process represented in Fig. [3]panel B), as well 
as the "re-combination" processes (see Fig. [3] panels C and D), are expressed in terms of the coefficients for the heavy-hole 
+3/2 reflection (panel A) given by Eq. llBl! . 



Finally, we define the symmetryzed eigenstates 

*l(r) = 
*k(r) = 
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W 2 



ii> BD+ ] 
iiP BD+ ] 



(19) 



where C\ and Ci are normalization constants and the linear combinations ip AC+ = ip A + ip c and i\) BD+ = i]: B + ip D 
contain wave-functions with coefficients given in Tabic IIIII The wave-functions defined by Eq. (|19J) represent a 
complete system of cigcnfunctions in the subspace characterized by k = |k| and the incident angle 9 > 9 C . The 
symmetry between heavy-hole contributions with opposite helicity is reflected by the special relation between the 
corresponding coefficients, C2 = ±i c\ and c§ = ±i C5. The complete set of relations between the coefficients is given 
in Table HV| 
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TABLE IV: Relations between the coefficients of the symmetric eigenstates in the presence of the localized modes, i.e. for an 
incident angle for the heavy-holes 8 > 9 C . The wave functions for the eigenstates fy 1 and given by Eq. 119t are expressed in 
the general form given by Eq. 118H . The coefficients b 2 and b\ are real. Notice that there are only two independent eigenstates 
for a given set of parameters (k, 8), in contrast with the case 8 < 8 C when four independent eigenvectors were found. 



We have identified complete sets of linearly independent eigenstates for arbitrary values of the incident momentum 
of the heavy-holes, i.e. for arbitrary parameters (k, 9). These sets of eigenstates are given by the equations <|14[) and 
(|19f) . In order to construct an orthonormal basis for the quantum problem defined by the Luttinger Hamiltonian in 
half-space, we have to solve two problems: a) the orthogonality of the eigenstates, and b) the proper counting of the 
modes. We address these issues in the remaining of this section. 



A. Orthogonality of the symmetrized eigenstates 

Let us start by clarifying a few aspects related to the definition of the scalar product that we use here. Consider two 
wave-functions and where k = (k x , 0, k z ), with k x /k z = — tan((9), represents the momentum of the incident 
heavy-hole mode and a E 1 — ,2+, 2—}, if 9 < 9 C or a £ {1,2}, if 9 > 9 C . For reasons that will become clear 
below when we address the problem of counting the states, we consider a "discretized" momentum space so that two 
momenta are considered identical, k = k', if they are both in a certain cell of volume 5k 3 and different otherwise. We 
define the scalar product of the two states as 



Jn 



(20) 
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where — L 3 is the volume of the system. We will always consider the thermodynamic limit L — > oo, so that L6k 3> 1. 
Due to the particular form of the eigenfunctions (|13[) and (|18fl . the integrations over x, y, and z can be performed 
separately. Due to our choice for the coordinate system, the wave-functions are independent of y and therefore the 
integration over this coordinate trivially generates a factor L. Further, all the x-dependence is contained in the factors 
e~xjp(ik x x) so that the x-integration generates a factor hSk x k' ■ Consequently, Eq. 1|20[1 reduces to the one dimensional 
integral 



(21) 



where ^f^(z) = ^^(r) exp(—ik x x). We notice that, in order to insure convergence, all the z-dependent oscillatory 
terms in Eq. (|21() are multiplied by a convergence factor exp(— rjz) with r]L 3> 1, and the limit 77 — s- is taken at the end 
of the calculation. Considering now the expressions (|13|) and Ijl8(l for the eigenfunctions, together with the properties 
of the coefficients summarized in Table ITT1 and Table Hvl as well as the expressions for the spinors given in Appendix 
1X1 we obtain the generic expressions of the terms contributing to the scalar product l|21|) . We will concentrate on the 
case k x = k' x7 because for different x-components of the momentum the states are manifestly orthogonal. Also notice 
that, due to the exponential factors exp(— Qz), the localized modes from the eigenfunctions (|18|l will generate terms 
proportional to 1/(LQ), which vanish in the thermodynamic limit, and therefore will not contribute to the scalar 
products. Omitting these terms, the generic contribution to the integral over z in the scalar product has the form 



L 2 [^l(z)]^i,(z)= Y, \ [Ar Q '(fci,fe)cos[(fc 1 -fc 2 )z]+A 2 l ' Q '(fc 1 ,fc 2 )sin[(fc 1 



(kiM) 



k 



2 Z 



(22) 



with (fci, fc 2 ) € {(k z ,±k' z ), (k z ,±q' z ), (q z ,±k' z ), (q z ,±q' z )}, where k z and k' z are the z-components of the incident wave- 
vector for the heavy-hole modes, while q z and q' z are the z-components of the incident light-hole waves. The reflected 
waves have components with opposite signs. Due to the symmetry properties of the eigenstates, some of the coefficients 
A"' Q are identically zero. The non-vanishing coefficients are summarized in Tabled Integrating now Eq. I|22|) over z, 





(1+) 


(1-) 


(2+) 


(2-) 


(1) 


(2) 


(1+) 


Ai 


A 2 


Aj 


A 2 


Ai, A 2 


Ai, A 2 


(1-) 


A 2 


Ax 


A 2 


Ai 


Ai, A 2 


Ai, A 2 


(2+) 


Ai 


A 2 


Ai 


A 2 


Ai, A 2 


Ai, A 2 


(2-) 


A 2 


Ax 


A 2 


Ai 


Ai, A 2 


Ai, A 2 


(1) 


Ai, A 2 


Ai, A 2 


Ai, A 2 


Ai, A 2 


Ai, A 2 





(2) 


Ai, A 2 


Ai, A 2 


Ai, A 2 


Ai, A 2 





Ai, A 2 



TABLE V: Non-vanishing A"' a (fci,fc 2 ) coefficients that contribute to the scalar product of the eigenstates , 3'k(- z ) an d ( z )- 
The rows of the table are indexed by a and the columns by a' . For example, the scalar product of 1 i> 2+ and will be an 
integral over z of a sum containing only |iamfoda 2 -type terms. The simple z-dependence of the integrand, given by Eq. J22J, 
and the properties of the coefficients A"' Q (fci, fc 2 ) enable us to prove the orthogonality of the eigenfunctions. 



and taking the necessary limit for the convergence factor 77, we obtain two types of contributions to the scalar product 
of two eigenstates ^ and 'J'^,: A"' a (fci, fci) if fci = fc 2 and j, A 2 '" (fci, k 2 )/(ki — k 2 ) if fci 7^ fc 2 . For the contribution 
proportional to A 2 we distinguish two regimes. The first regime is characterized by a finite difference fci — fc 2 , so that 
L\k\ — k 2 \ 3> 1, and occurs for example in combinations like k z + k z or k z + k' z . In this case there is no contribution 
to the scalar product in the thermodynamic limit. In the second regime, the difference between the two components 
of the wave-vector can be arbitrarily small, so that fci ks fc 2 and the ratio A 2 ' a (fci,fc 2 )/(fci — fc 2 ) seems to diverge. 

However, by evaluating explicitly the coefficients in this limit we obtain A^'™ (fci, fc 2 ) cx fci — fc 2 and the corresponding 
contribution to the scalar product vanishes again in the thermodynamic limit. This behavior is valid for all cases, 
except when (a, a') = (i±, j=p) with i 7^ j, for example in the case k z w k' z , q z « q' z and (a, a') = (1+, 2—). In these 
special cases when the A 2 terms do not vanish individually, a direct evaluation of the coefficients in the limit k z — > k' z 
shows that 



A 2 '° (q z ,q' z ) 
A 2 ' (k z ,k' z ) 



dk z 



C|cos(g)| 
C- sin 2 (9) 



(23) 
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Consequently, the total contribution to the scalar product, -^[A^'" (k z , k' z )/(k z — k' z ) + A^'" (q z , q' z )/(q z — l' q )] vanishes 

again in the thermodynamic limit. The A"' a (fci,fcx) terms also vanish for a ^ a'. In fact, the only contributions to 
the scalar products that survive in the thermodynamic limit arise from A"' Q (fc 2 , k z ) and A"' a (g z , q z ) and we have 

8 

A? a (k z ,k z ) + A«- a (q z ,q z )=Y, |cf | 2 = 1 for 9 < 8 C , 

i=l 

h-TikM = K\ 2 + \c%\ 2 + \ct\ 2 + \c%\ 2 = 1 for 0>0 C , (24) 

where the unity was obtained by a proper choice of the normalization constants for the coefficients cf in the con- 
struction of the symmetrized states. We conclude that the symmetrized eigenstates constructed here represent an 
orthonormal system satisfying (^kl^k') = ^aa'4k'- 



B. Counting the eigenstates 



The standard procedure for counting the states and transforming sums over the k- vector into integrals for a quantum 
system in the thermodynamic limit consists in imposing periodic boundary conditions for a finite volume f2 = L 3 
system, and then taking the limit L — > oo. The periodic boundary conditions generate a discrete set of wave- vectors 
k = (k x , k y , k z ) = 27r/L(ni, ri2, 713), where n, are integers. The procedure works in the presence of spin-orbit coupling 
if no boundary is considered and the eigenstates arc the heavy-hole and light-hole modes. However in our case, the 
eigenstates are certain combinations of heavy- and light-holes corresponding to the same energy. Explicitly, for 9 < 6 C , 
heavy- hole modes with k# = 2it / L(n x , n y ,n z ) have to be paired up with light-holes characterized by the wave- vector 
(\l = 2tt / L(n x ,n y ,m z ) with m z satisfying 

f (t) 2 {nl + n y + *«) = (t) 2 {nl + n " + ml) ■ (25) 

For an arbitrary value of £ = mi jvriH there is no integer solution of this equation. To overcome this difficulty we 
consider a partition of the reciprocal space with the property that all the heavy-hole modes characterized by helicity 
A and wave- vectors situated in a cell Sk 3 centered on k are represented by a single heavy-hole mode (k, A) with the 
energy = k 2 /2ttih and an effective "degeneracy" i/q = <5fc 3 (L/27r) 3 ^> 1. In the thermodynamic limit we can choose 
8k arbitrarily small, so that 1/Sk is larger than any relevant length-scale in the problem. For incident angles 8 < 6 C , 
the heavy-hole modes from a cell (k, 5k 3 ) will pair with light-hole modes from a cell (q, Sq 3 ), were all the points of the 
new cell were obtained by mapping the points of the original cell, k — > q, using Eq. (|10f> . The ratio of the light-hole 
and heavy-hole degrees of freedom from the two cells is 



8q 3 dq z t\cos{6)\ 



Sk 3 



dk. 



\A-sin 2 (0) 



(26) 



For each heavy- hole mode with helicity +3/2 used in the construction of the eigenstates 1)14(1 we will use another 
heavy- hole with opposite helicity as well as 2Rlh light-hole modes and we will generate 2/xi states ^E 1 ^ and 2/X2 
states V I' 2± . The effective degeneracies ^ are determined by the condition that, far from the boundary, the number 
of heavy- and light-hole modes is the same as in an infinite system. Consequently we have 

rtl + |tt^ 8 = 1, 

HiWl% + i*Wl% = Rlh, (27) 

where Wlh and Whh are the light-hole and heavy-hole weights, respectively. We obtain for the effective degeneracies 
of the eigenstates ^Jf given by l)14|l the expressions 

Ml ~ WFh - W¥n ' (28) 

Notice that these expressions are meaningful only if the effective degeneracies satisfy the condition < /Xj < (1 + Rlh)- 
Consequently, the weights of the heavy-hole modes contained in eigenstates defined by equations l|14|) cannot be 
arbitrary. If, for example, W][ H > Wj[ H they have to satisfy the inequalities 

Whh > ~ , W&b < . (29) 
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These inequalities are the reason behind our choice for the parameter 'a' in the construction of the eigenstates i|14f> . 
By maximizing the weight of the heavy- holes in the type-1 eigenstates and minimizing it in the type-2 eigenstates we 
insure that the inequalities (|29J) are always satisfied. The analysis of the case 9 > 9 C is straightforward, because far 
from the boundary only the heavy-hole modes survive. Consequently, there is a direct correspondence between the 
number of degrees of freedom associated with the heavy- holes and the number of symmetrized eigenstates <|19|) . 

In this section we have constructed an orthonormal basis for the Luttinger Hamiltonian in half space. For our 
particular choice of coordinates, the basis consists in the symmetrized eigenfunctions given by equations i|14|) 
with i = 1,2, <r = ±, and k = (k x ,0,k z ) for incident angles 9 — — arccos(fc 2 /fc) < 9 C , together with the eigenstates 
given by Eq. i|19[l with i = 1, 2 and incident angles 9 = — arccos(fc z /fc) > 9 C . The eigenstates from the first set 
have an effective degeneracy fo/^i with /i^ given by (|28|l . while the eigenstates from the second set have an effective 
degeneracy Vq in our discretization construction for the momentum space. Notice that the overall factor v is irrelevant 
in the calculation of the physical quantities. We conclude this section by giving the rules for transforming sums over 
the k-vector into integrals. Let us assume that f(k, a) is a function that depends on the diagonal matrix elements 
(^gm^g) of a certain operator A and that we are interested in calculating its sum over all the eigenstates of the 
basis. We have 

£ 5>0Mi)/(k,*» + |5 Hm)/(k,0 = ^X; / d 3 k l i i f(k,i,a) + -^Yl I d 3 kf(k,i), 
k it k i (^)~J(e<e c ) N iW) 

(30) 

where we have taken into account that the volume of a cell in the discretized momentum space is Sk 3 = vq(2it/L) 3 
and we have considered the thermodynamic limit L — * oo. For example if /(k, a) represents the occupation number 
at zero temperature of the mode with the quantum numbers (k, a) , by applying Eq. (|30fl we can easily obtain the 
relation between the Fermi k-vector and the density of particles, n = k 3 F {\ + ^i)/(6Tr 2 ). Notice that each component 
^ of the basis, although contains a superposition of heavy- and light-holes, is labeled by the k-vector of the incident 
heavy- hole mode and has an energy = fc 2 /(2m#) independent of a. Consequently, the Fermi wave- vector kp is 
independent of a and is determined by the density of particles and the spin-orbit coupling. 



III. OSCILLATING SPIN DENSITY AND SPIN ACCUMULATION 



In this section we investigate the properties of the spin density and explore the possibility of spin accumulation in a 
system with a sharp boundary described by the Luttinger model (here and below we actually study the accumulation 
of the total orbital momentum, but we call it "spin accumulation" for brevity). We can derive the spin density from 
the one-particle propagators. Using the orthonormal basis constructed in the previous section, we can define the 
Green's functions for the Luttinger model in half space. The retarded Green's function is 

G ( *W r rW W mrTb in{r)l (31) 

(w ' r ' r )-2^ w _ /x _ efe + z0+ < ( 31 ) 

k,a 

where /j, is the chemical potential and ^(r) ^ s a vec t° r from the basis defined by the equations Ijl4[l and i|18[) 
corresponding to the eigenvalue = k 2 /2itih, with run being the heavy-hole mass. The eigenvectors are labeled by 
the wave-vector k of the incident heavy-hole mode and the index a which takes four values, a g {1+, 1 — , 2+, 2—}, 
if all the modes propagate in the z-direction (9 < 9 C ), and two values, a € {1, 2}, if the cigenstate contains localized 
modes (9 > 9 C ). In addition, because the states 4 , £( r ) are four-component spinors, we have the indices a and b to 
label their components. Finally, to correctly account for the light- and heavy-hole modes used in the construction of 
the basis, we have the effective degeneracy factor v^ a = VQii a , where \x a is given by Eq. I|28() when the eigenvector 
contains only propagating modes, i.e. for 9 < 9 C , and /i a = 1 otherwise, while isq = Sk 3 / (2ir / L) 3 is determined by our 
momentum discretization procedure (as we have noted above this factor drops out of all final results). The density of 
the total momentum S can be expressed in terms of Green's function as (S)(r) = —i J du)Tr[SG(u; r, r)], where the 
matrices S a f, = ([S^lab, [S y ] a b, [S z ]ab) are given in Appendix lAl In terms of the eigenstates for the Luttinger model in 
half-space we have explicitly 

(«)(r)=EE^ ([*kto] + «a(k), (32) 

k a 

where Si is the matrix for the I € {x, y, z} component of the total momentum, n Q (k) € {0, 1} is the zero temperature 
occupation number of the ^ state with effective degeneracy z^ Q , and the summations, which include all the eigenstates 
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in the basis, can be converted into integrals using Eq. (|30f) . The analysis can be greatly simplified if we take into 
account the symmetry properties of the eigenstates. These symmetries translate directly into relations between single 
state contributions to the matrix clement of the spin density, 

5.r(k,r) = [*£(r)] t 5,vl/£(r). (33) 

Explicitly we have 

|>Sr(k,r) = -5i M (k,r) for 9 < C , 

\ SI (k,r) = for 9 > 6 C , 

f<Si CT (k,r) = for 0<9 C , 

\ S l z (k,r) = for 9 > C , 

where i = 1, 2 and a = ± represent the labels of the corresponding eigenstates. We conclude that the contribution 
to the x and z-components of the spin density from any wave-vector k is identically zero. On the other hand, for the 
y-component we obtain 

for 6 < 6 C , 
for 9>6 C , 

and consequently we get in general a non-zero total contribution for a given k. In other words, in a system described 
by the Luttinger Hamiltonian and having a sharp planar boundary, the contribution to the spin density vector from 
states characterized by a wave- vector k is oriented along a direction parallel to the boundary and perpendicular to 
the wave-vector, <S(k, r) oc nj x k, where n& is a unit vector normal to the boundary. We remind here that the 
cigenfunctions ^ arc written in the "local" coordinate system determined by the convention k = (k x > 0, 0, k z ), and 
therefore, whenever we have to use a system of coordinates independent of k, special attention should be given to the 
transformation rules. In particular, for the spin density vector we have 

S((k x ,k y ,k z ), r) = -S((-k x ,-k y ,k z ), r), (36) 

i.e. the contributions to the matrix element l|33|l from wave-vectors with opposite components parallel to the surface 
have the same length and opposite orientations. Consequently, in the ground-state the spin-density vanishes. However, 
it is possible to generate a non-vanishing spin-density if the number of particles moving in one direction is different 
from the number of particles moving in the opposite direction, i.e. in the presence of a charge current parallel to the 
boundary. 

Let us assume now that the system is characterized by a preferential direction of motion for the carriers imposed 
by an external current flowing in the x-direction and having the average current density (j x )- We reiterate that if the 
system is perfectly clean and an external electric field is present, then the particles will just accelerate indefinitely 
and the spin density will be timc-dependent. This scenario is not considered here. Below we study the situations 
when the equilibrium current appears either due to a voltage drop which occurs entirely in the contacts, or because 
some disorder is present. In the latter case, we study only the lengthscales much smaller than the mean-free path (the 
physics at larger length-scales in disordered systems will be discussed in Sec. IIV|I . We can associate the current with 
a drift velocity vo = (Ak/rriH, 0, 0) with Ak <C fc^, and assume that the distribution function for the system changes 
from n°(k) = Q(cf — £k) to n(k) = 0(ej? — £k+Ak), where 0(x) is the step function, et = k 2 /2m# are the energies 
of the eigenstates, and €f is the Fermi energy corresponding to the Fermi wave- vector kp related to the density of 
particles, n ~ k F (l + £2 )/(67r 2 ). Using the relation l|3()(l of transforming the sums over k into integrals, we have 

0*> = A— I / d " k t 1 + Rlb(*)] fc*n(k) + / d 3 k k x n(k)\ = nAfc, (37) 

47T J m ff lJ e< g c J 8> g c J m H 1 + £2 

where n is the average particle density of the system and we have taken into account that the drift velocity is much 
smaller than the Fermi velocity. In the presence of the external current, the y-component of the spin density becomes 
non-zero and, using the definition l|32|) together with Eq. (|30ll . we have 

WW - J^y ^ ^msl a (O,r)+V2(O)S 2 y °(0,r)] cos(0) »(k) (38) 

r ~) k 2 Ak { f 6c r /2 1 

+ J d 3 k [5i(0,r)+5 y 2 (0,r)]cos(^n k j-^-^M o d6S < {e,z) + J^ d0<S>(0,z)L 



13 




where Hi is the effective degeneracy of the eigenstate, ra(k) = 8(ef — ek+Ak), <f> is the angle between the x-axis and 
the projection of the wave- vector on the x-y plane, and the factor cos(0) is due to the fact that S y (9, r) is calculated 
in the local coordinate system. Finally, the functions integrated over the angle 9 in the last member of Eq. Ij38(l are 

S<(0, z) = 2 [m{0)Sl + (9, z)+^ 2 (9)S 2+ (9, z)} sin 2 (6) and S > (9, z) = [S^{9, z) + S 2 {9, z)] sin 2 (9). (39) 

As expected, the spin density is proportional to Ak, i.e. to the average density of the external current. So far, we 
have reduced the expression for the spin density (S y )(r) to a sum of one dimensional integrals over the incident angle 
9. Furthermore, we can write explicitly the z-dependence of S(9, z) using the expressions for the wave-functions and 
spinors as well as the symmetry properties of the coefficients. We obtain 

S < (9,z) = B{6) ism[2k z z] + $m[2q z z] - 2sm[(k z + q z )z]\ , 

S>{6,z) = B{9) |sin[2fc 2 z] - 2 sin^e^ 2 } +C{9) \cos[2k z z] + e -2 ^ - 2cos[fc*z]e -(5 *} , (40) 

where the wave- vectors are proportional to hp and depend on the incident angle 9. Explicitly we have 

k z = fci?cos(6*), 

q z = k F ^-sm 2 (9), 9<9 C , (41) 

Q = k F ^sm 2 {9) 9 > 9 C . 

The complete information about the spin density is contained in the coefficients B(9) and C(9). These quantities can 
be expressed in terms of the coefficients cf of the symmetric eigenfunctions (|14|) and (|18|l . Explicitly, we have for 
9<9 C : 

B(9) = -6 {/xi [(c\) 2 - (cl) 2 } + M2 [(c 2 ) 2 - (c 2 ) 2 ] } cos(0) sm 4 (9) - 12 [^c\c\ + ^cl] cos 2 {9) sin 3 (0), (42) 

where fj-i are the effective degeneracies given by equation (|28|l and the notations for the coefficients cf are taken from 
Table ITT1 Similarly, for 9 > 9 C we have 

B{9) = -6Rc [(c}) 2 + (c 2 ) 2 ] cos(6») sin 4 (6>) + 6Rc [i{c\f - i(c 2 ) 2 ] cos 2 (6») sin 3 (6i) 

C{9) = 6Im[(c}) 2 + (c 2 ) 2 ]cos(0)sin 4 (0)-6Im[z(cJ) 2 -i(c 2 ) 2 ]cos 2 (0)sin 3 (0), (43) 

where the notations for the coefficients cf are taken from Table IfVl The explicit dependence of the coefficients B and 
C on 9 is shown in Fig.[3]for several values of the mass ratio £. Notice that the coefficient C(9) vanishes at the critical 
angle 9 = 9 C and at ir/2, while B(9) vanishes at 9 — 0, 7r/2 and has a singularity at the critical angle characterized 
by a discontinuity of the first derivative. Knowing these coefficients, we can perform the integral over the angle 9 in 
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FIG. 6: (Color online) Position dependence of the current-induced spin density (S y ) expressed in units of k'pAk /8n 2 for a mass 
ratio £ = 0.3 relevant for a GaAs system. The inset shows the large-z behavior characterized by a decay of the oscillations 
proportional to 1/z 2 (blue lines). 



Eq. (|38|) and obtain the spin density. The result for £ = 0.3 is shown in Fig. revealing several interesting features. 
It confirms that the current-induced spin density is indeed non-zero. In fact the spin density oscillates about zero 
with an amplitude that vanishes as we move away from the boundary. The oscillations are characterized by several 
wave-lengths, as shown by the beat phenomenon in the inset of Fig. [5] On the other hand, the amplitude of the 
oscillations decreases as 1/z 2 at large distances. All these features are contained in the properties of the coefficients 
B(9) and C(9), in particular in their behavior in the vicinity of the singular points 9 = 0, 6 C and tt/2. Below, we will 
analyze in detail the connection between these singular points and the oscillating properties of the spin density. 



where 



A. The asymptotic behavior of the spin density 

The amplitude and the period of the spin density oscillations far from the boundary are completely determined by 
the singular points of the coefficients B(9) and C{9). To prove this, let us consider as an example one of the terms 
from Eq. I|40|l and re-write the corresponding contribution to the spin density as a Fourier transform, 

e a poo 

dO B(6) sin(2k z z) = / dp b (p) sin (pz), (44) 

J — OO 

!B[ arc cost , p ) ) , 
-i ; {2kF,) if 0< P <2y/T=T, 
~ (45) 
otherwise. 

If the Fourier coefficient b(p) has a singularity at po ^ characterized by a discontinuity Ab^(po) i n the derivative 
of order n, the large-z behavior of the integral l|44Jl is Ab^ n \po) cos(poz + n-n /2) / z n+1 , i.e. it oscillates with a period 
A = 2ir/po and an amplitude which is proportional to the discontinuity in the derivative and decays as l/z n+1 . Similar 
relations can be established for cases when the discontinuity in the derivative is not finite. For example, when the 
singularity in proportional to \fp — po, the asymptotic behavior becomes cos(poz + 7r/4)/z 3 / 2 . Returning now to the 
analysis of the contributions to the spin density coming from Eq. I|4(J|) , we notice that the possible singularities are 
given by the angles 9 6 {7r/2, 9 c , 0}. As k z (9 = tt/2) = 0, there is no large-z contribution to the spin density coming 
from the states with the incident angle 9 = tt/2. In other words, the modes propagating parallel to the boundary do 
not contribute to the spin density at positions far from the surface of the system. On the other hand, the critical angle 
9 C requires a detailed evaluation. We should remind that, as shown in the inset of Fig. El the spin density oscillation 
decay as 1/z 2 . From Eq. (|40|) we observe that there are eight different terms that contribute to the spin density, 
three coming from 5< and five coming from 5>. However, one can show that the contributions coming from the last 
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FIG. 7: (Color online) Total spin accumulation in units of kpAk/8n as a function of the mass ratio £ (black line with triangles). 
The contributions from angles 8 < C (orange line with diamonds) and 6 > 9 C (green line with circles) are shown separately. 
The inset represents a fit of the total spin accumulation (triangles) with the function (1 — £) 2 (continuous red line). 



two terms proportional to C(9) decay faster that 1/z 2 and, consequently, can be neglected. The remaining six terms 
will give contributions that decay as 1/z 2 or slower. For example, the first term B(9) sin(fc 2 z) with 9 < 9 C generates 
a contribution proportional to B(9 C ) cos(2^/l — £z)/z due to the fact that B(0 C ) ^ 0. However, this contribution is 
exactly canceled by the corresponding contribution coming from the fourth term, B(9) sm(k z z) with 8 > 9 C . Similarly, 
the a/6*c — 9 term from B{9) generate the contribution cos(2\/l — £z + 7r/4)/z 3 / 2 , which is exactly canceled by the 
contribution generated by C{9). In fact, one can show that there is no net contribution to the asymptotic spin density 
of order up to 0(1/ z 2 ) coming from the critical angle. This exact cancellation is embedded in the properties of the 
coefficients B(9) and C(9) in the vicinity of 9 C . Explicitly, for 9 sw 9 C we have 

B(9) = 0o + p 1/2y /0 c -0 + l3 l (9 c -8) + O ((0 C - 9) 3 ^ , 9 < 9 C 

B(9) = fa - Pi(0 - 9 C ) + O ((8 - 9 C ) 2 ) , 8>8 C (46) 

C{9) = 1/2 V6-Oc + O ((#-0c) 3/2 ) , 0>0 C . 

We conclude that the asymptotic behavior of the spin density is determined entirely by the modes propagating along 
the direction perpendicular to the boundary. The characteristic wave- vectors corresponding to the first three terms 
in Eq. igUJ) are 2fc z (0) = 2k F , 2q z (0) = 2k Fy /£ and L(0) + q z (Q) = k F (l + V?)- On the other hand, the amplitude is 
determined by the small angle behavior of the coefficient B, namely 

B{9) =-39 3 + 0(9 5 ), for 9^0. (47) 

Writing each term as a Fourier transform and using the properties discussed above in connection with equations 
()44I45JI we obtain the following expression for the spin density far from the boundary 

(S v )(z) = ^ ^ | sin (2k F z) + gsin (2k F y/lz) - sin [k F (l + ^)z] j , zk F > 1, (48) 

where the Fermi wave-vector is determined by the density and the spin-orbit coupling, fcp = [6n 2 n/(l +^ 3 / 2 )] 1 / 3 . We 
note that the result l|48() is asymptotically exact. 



B. The spin accumulation 

An important question is whether there is a net spin accumulation in the system. To answer this question, we first 
define the spin accumulation as an integral of the spin density with respect to the distance from the boundary 



Sy(z) = / (S v )(0 dC (49) 
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FIG. 8: (Color online) Current-induced spin accumulation in units of kFAk/8n as a function of the distance from the boundary. 
The curve was determined for £ = 0.3 by integrating the data from Fig. HJwith respect to z from zero to a certain distance from 
the boundary. Notice that the spins are localized within a few Fermi lengths Af = 2n/kF off the boundary. The dashed line 
represents the total spin accumulation S y using the same units. The inset shows the 1/z 2 decay of the fluctuations about S y 
at distances far from the boundary. 



Here S y (z) is the spin accumulation at a distance z from the boundary, more precisely the spin density per unit area. 
The total spin accumulation is S y — S v {oo). To determine the total spin accumulation it is convenient to perform 
first the integration over the coordinate z for each of the terms contributing to S(9, z). We obtain for the total spin 
accumulation the expression 



S n 



d9 B{9) 



1 

2kZ 



1 

Wz 



k z + q z 



tt/2 



d6 {B(6) 



1 

2kZ 



k 2 z + Q 2 



+ C(0) 



1 

2Q 



2Q 



(50) 



where the dependence on 9 of the wave- vectors k z , q z and Q is given by Eq. I|41|l . Finally, the integration over the 
angle is performed numerically. The result for the spin accumulation as a function of the mass ratio £ = rnj,/mir 
is shown in Fig. {7\ The spin accumulation decreases monotonically with £ and vanishes in the absence of spin-orbit 
coupling (£ = 1). We have also calculated separately the contributions from the angles smaller than the critical angle 
9 C (orange line in Fig.Ql and the angles larger than 9 C (green line). For strong to intermediate values of the spin-orbit 
interaction practically the entire contribution to the spin accumulation comes from states with 9 > 9 C , i.e. the states 
containing localized modes. In contrast, in the weak-coupling regime (£ ~ 1), the two contributions are comparable 
and have opposite signs. Finally, we notice that the ^-dependence of the total spin accumulation can be perfectly 
fitted as (see the inset in Fig. 0) 



(51) 



The spin accumulation as a function of distance has also been calculated numerically and the result is shown in 
Fig. [H] We see that the spin accumulation is localized in the vicinity of the boundary. The shortest distance from the 
boundary at which the spin accumulation becomes equal to the total spin accumulation S y is of the order of 2/fcp, 
while at larger distances it oscillates about S y with an amplitude that decreases as we move away from the boundary. 
As shown in the inset of Fig. [SJ at large distances the oscillations of the spin accumulation decay as 1/z 2 . 



IV. THE EFFECTS OF DISORDER 



A. The Friedel oscillations in a disordered system 



In this section we qualitatively discuss the effects of disorder on the current-induced Friedel oscillations of the spin 
density. To understand these effects it is useful to recall the corresponding physics in the usual disordered Fermi 
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systems: Consider a perturbation (such as a boundary or an impurity) introduced in a Fermi gas. This perturbation 
leads to the Friedel oscillations in the electron density, which in a clean case behaves as n (r) = Ar-P^ cos (2p F r), 
where r is the distance from the perturbation, A is a constant amplitude, and the power p(d) = (d + l)/2, if the 
perturbation is a boundary and p(d) = d, if the perturbation is point-like (i.e., a magnetic impurity). What happens 
with the Friedel oscillations if we introduce disorder? A natural quantity to calculate is the perturbation-induced 
density averaged over disorder realizations. As was shown by de Gennes^ it decays exponentially as (n(r)) oc e~ r l l , 
were / is the mean free path. However, Zyuzin and Spivak^a later argued that this result is misleading. Indeed, if 
we measure the density at different points but at the same distance away from the perturbation, we will get different 
results (in magnitude and sign), which will strongly depend on the local disorder realization (see also, Refs. l26j| ) . 
Therefore a more meaningful quantity is the full probability distribution function of density at a distance r from the 
perturbation P[n,r]. One should note that this probability distribution function is not known explicitly even in the 
simplest cases and only its lowest moments are available. The important result of Zyuzin and Spivak24 (obtained in 
the context of RKKY oscillations between magnetic impurities in disordered metals) concerns the second moment, 



which was shown to decay as a power law, (n(r)n(r)) 



-P(d) \ 



just like in the clean case. The qualitative explanation 



of this result is that disorder effectively introduces random phase shifts 5(1) in the argument of the cosine in the 
Friedel oscillations cx cos [2p-pr + 5(1)]. Averaging over disorder realizations leads effectively to the vanishing density 
due to the oscillations in the cosine. However, averaging of the cosine squared yields a finite result and restores the 
power law behavior of the typical Friedel oscillation. 

To connect the above-mentioned arguments with our case, we mention that the higher moments of the probability 
distribution function of the Friedel oscillations in the density are described diagrammatically by bubble diagrams 
with all possible disorder averagings between them. A careful treatment of these diagrams reveals that the power-law 
(non-exponential) decay of the usual Friedel oscillations is due to the gapless nature of the diffusons of which these 
diagrams are built. Now let us return to the Friedel oscillations in the spin density and spin accumulation, which were 
discussed above in the clean case. A straightforward calculation shows that the spin density averaged over disorder 
decays exponentially as 



S(r) oc / deTr SG(e;r,r 



cx e 



-r/t 



(52) 



where I is the mean free path, G(e; r, r) is the Green's function which takes into account the boundary, and overlinc 
implies disorder averaging. Higher moments of the spin density are described by expressions, which involve different 
components of the diffuson. In particular, the second moment (correlator of spin densities) reads 



S Q (r)S Q (r) oc J de x J de 2 Tr [s a G( £l ; r, r)] Tr [s Q G(e i; r,r) 



1 



(53) 



We emphasize that due to the two traces present in Eq. (|53|l . the correlation function includes contributions from 
different components of the diffuson matrix. Most importantly it contains a singlet component, which remains gapless 
even in a spin-orbit coupled system (see below subsection II V Hj) . Therefore, the typical spin density measured at a 
distance r I from the boundary is expected to decay as a power law even if disorder is present. This spin density is 
however random in sign and can be estimated as S y (r) oc z~ 2 sin [0 r (r)], where z is the distance from the boundary 
and cj> r is a random phase. Due to this randomness, the spin accumulation (i.e. the spin density averaged over the 
length-scales much larger than the Fermi-wavelength) is expected to decay exponentially as oc e~ z ^ Ls , where L s is the 
spin relaxation or spin diffusion length>i2i2L2& We note that in the framework of the Luttinger model, the spin-orbit 
coupling is always strong and therefore the spin relaxation length is expected to be very short and of the order of the 
mean free path. We study the related spin relaxation time in the following subsection. 



B. Spin relaxation time in the Luttinger model 

For the sake of completeness, we calculate explicitly the spin relaxation time and find that it is very short (of order 
r). Therefore, the spin relaxation length is short too and the hydrodynamic diffusion approximation is not applicable 
for spin transport in the Luttinger model. However, we believe that it is likely that the result for the spin relaxation 
time itself is quantitatively correct and describes time relaxation of a uniform spin distribution. 

To find the spin relaxation time, we construct the diffuson, which can be obtained from a convolution of disorder- 
averaged Green's functions [see Eqs. (|57|l and l|58|l below]. For the purpose of this section, we can ignore the boundary 
effects and use the bulk Green's function in the calculations. We consider the system in the presence of randomly 
distributed impurities with short-range potential Vi mp (r) = uq 5(r — r), and calculate the self-energy in the Born 
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approximations*. We obtain for the retarded (advanced) self-energy the expression 



^kA,k 



v(w±0+) 



2r 



(54) 



where pp — miikpil + £, 3 ^ 2 )/tt 2 is the density of states at the Fermi energy, m represents the density of impurities, 
and r = 1/ ' {jmiV^pp) is the mean scattering time. Notice that the self-energy in the first Born approximation is 
momentum independent and diagonal in both momentum and helicity. We consider here the quasi-classical limit, 
when the mean- free paths for heavy- and light- holes are much longer than the corresponding Fermi lengths, l\ = 



k^'r/mx ^> 2Tt/kp . Using the expression (|54|) for the self-energy, we obtain for the disorder-averaged retarded 
Green's function the usual expression^ 



f R (k,w) + ^g R (k,u,) 



/ - g R (k, w)(n k S) 2 = f R (k, lo)I + g R {k, to) d(n k )I\ 



(55) 



where n k = k/|k| is the unit vector directed along the momentum, / represents the 4x4 unit matrix, the S-matrices 
are given in appendix lAl and d(n k )r = dj(n k )l\ represents the product of a five-component unit vector that depends 
on the direction of k. In what follows we use the technique introduced by Zhang et ai, which involves the five 
generators l\ of the SO (5) Clifford algebra^. Finally, the coefficients f R and g R are independent of the direction of 
the momentum and read 
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(56) 



where p is the chemical potential. Similar expressions can be obtained for the advanced Green's function with 
f A {k,w) = [f R {k,w)Y and g A (k 1 uj) = [g R (k,u)}*. 

We are ready now to calculate the kernel of the diffusion equation, 



v = [/® j-nj- 



(57) 



where T> and n are 16 x 16 matrices labeled by two pairs of "spin" indices a <E {3/2, 1/2, — 1/2, — 3/2}, and I is 
the 4x4 unit matrix. The polarizability can be expressed as a convolution of two Green's functions, 
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d 3 k 



■nppT J (2tt) 
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(58) 



Considering now the limit lo <C ep and q = in Eq. (|58(l and performing the integral over momenta we obtain 



I%,0) = -(1 



IUJT) 



I® I- 



+ ^K)(1 + »wt) 



I® I- 



(59) 



where / is the 4x4 unit matrix, fj are the gamma- matrices defined in |T^ . and the coefficient T depends on the 
strength of the spin-orbit coupling. Explicitly we have 



^(0 



4(epT) 2 l + e/2 (l_^ + a ^(l+0 2 



(60) 



where £ = m^/m^r = (1 — 27)/(l + 27) is the light-hole to heavy-hole mass ratio and ep is the Fermi energy. Notice 
that epT ^> 1 and, consequently, the coefficient T is small, T ~ 0(l/(epT) 2 ), for any value of the spin-orbit coupling 
satisfying the inequality 7 ^> 1/(4eft), which is always the case for physically relevant parameters. We recall that 
the small spin-orbit coupling limit is not physically consistent within the Luttinger model (if 7 — > 0, additional bands 
have to be considered), although mathematically well-defined. Formally, we have .F(l) = 1/2 and we recover the 
standard result for systems without spin-orbit coupling. If we substitute now the expression (|59|l of the polarizability 
into Eq. JS7J, we obtain for the diffuson 



v- 1 = 



1® 1 



r, > — iujr 



i®i+ 



(61) 
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The final step in our evaluation is to write the diffuson in a different representation that has a more transparent 
physical meaning. Explicitly, we define 



V a b — M% i<72 V aiC j 2trT3aA M. 



b 

0-30-41 



(62) 



where M% a are 4x4 matrices and summation over repeating indices is implied. The M-matrices satisfy the resolution 



of identity M« 1(T2 M£ ( 



an d have the property Tr[(M a ) 2 ] = 1. In particular, 



j/2 represents 



the charge channel, while M 1 = S x /y/Z, M 2 = S y ^/5, and M 3 = S z ^/5 correspond to the spin channels. Notice that 
the original quantity is a 16 x 16 matrix and, consequently, there will be 12 other channels, in addition to charge and 
spin, corresponding to higher magnetic moments^ In fact Eq. (|62J) can be viewed as a multipole expansion of the 
density matrb*^ with the additional channels associated with the five T-matrices which are quadratic combination 
of the spin matrices^ ( "quadrupole moments" ) and seven linearly independent matrices containing cubic termsi^ 
("octupole moments"). Introducing the result (|61|l in Eq. (|62() we obtain the kernel of the diffusion equation in the 
new representation 



^>,0) = 
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(63) 



with 



a a = 1 - - 
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Tr[M a t i M a T. l 



Pa = l + -^Tt[M°f i M a f i 



(64) 



Finally, from equation (|63|l we evaluate the relaxation times 

[±-^(0K 1 
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a a 



±/3 a +^(£K T p a r' 



(65) 



where for the approximation we considered physically relevant parameters satisfying 1 — £ 3> l/(e_pT). Formally, in 
the limit of vanishing spin-orbit interaction £ = 1 we have — 1/2 and all the relaxation times become infinite. 
However, in the physically relevant regime, relaxation is naturally absent only in the charge channel as ag = and 
Po = 2. The spin relaxation times are 



The relaxation times for the "quadrupole" and "octupole" channels are 

3 



and T n 



Notice that all these relaxation times do not depend on the spin-orbit coupling [if 1 — m^/m^r S> (eft) 1 ] and arc 
comparable with the scattering time. Thus, these channels can not be described within the diffusion approximation. 



V. CONCLUSION 



We have found that the boundary spin Hall effect in a clean hole-doped semiconductor can be viewed as current- 
induced Friedel oscillations of the spin density with a net spin accumulation originating, in particular, from the 
localized surface states of the light-holes. Although, the explicit results derived in this paper are specific to the 
Luttinger model, the main qualitative conclusions generalize to other spin-orbit coupled systems^*^ Indeed, the 
existence of a Fermi surface in a metal or a semiconductor inevitably leads to the Friedel oscillations if a perturbation 
(such as a boundary or an impurity) is introduced in the Fermi gas. In the presence of a spin-orbit coupling, there 
are multiple Fermi surfaces and therefore there should be Friedel oscillations with multiple periods. It is also clear 
that the existence of multiple Fermi surfaces should lead to the appearance of surface states: As one can see from the 
derivation of Sec. ^ these states occur due to the impossibility to satisfy the energy and momentum conservation 
laws simultaneously if a majority carrier is reflected from the boundary having a large angle of incidence. This results 
in the appearance of localized surface states of minority carriers^*^ We note that since different types of carriers 
appear due to the spin-orbit splitting and have different spin properties, the surface states of minority carriers will 
always lead to a boundary spin accumulation if there is an equilibrium charge current present. These boundary 
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states may also be important in the context of spin diffusion in disordered systems with weak spin-orbit coupling. 
Indeed, it has been shown that the boundary spin accumulation in such systems is extremely sensitive to the boundary 
conditions i^2Si2L22i2£ The existence of the surface states may affect boundary conditions, since the "localization" of 
the minority carriers at the boundary would mimic a partially transparent interface. We note here that this effect 
is definitely non-perturbative and cannot be captured within expansions in powers of spin-orbit coupling, which is 
usually used in the context of spin diffusion. 

It seems plausible that the existence of the Tamm-like surface states can be connected to the Fermi surface Berry's 
phase structure in spin-orbit coupled systems, but this connection remains unclear at this stage. It is also unclear 
whether the spin Hall conductivity can be related to any of the observable properties in the spin Hall type experiment. 
In our analysis, we have not been able to identify such quantities, although the orientation of the spin accumulated 
near the boundary qualitatively agrees with the predictions of the spin-current theory. The main inconsistency is that 
the spin Hall conductivity tends to a finite constant if the spin-orbit coupling is small, while the spin accumulation 
calculated here depends strongly on the spin-orbit coupling 7 and vanishes if the latter goes to zero. 

We emphasize that our qualitative analysis is applicable only at length-scales much smaller than the mean-free 
path. At larger length scales different mechanisms of spin accumulation come into play. One such mechanism is due 
to the Dyakonov-Perel-like spin relaxation^Siii^ which together with the spin-charge coupling may provide another 
contribution to spin accumulation. However, as shown in Sec. II VI the spin relaxation times are always too short in 
the framework of the Luttinger model and therefore the above-mentioned disorder-induced spin accumulation will 
decay exponentially at distances of order mean free-path. In particular, this implies that the semiclassical diffusion 
approach is not applicable to the problem. 

In summary, we have constructed an exact orthonormal basis for the isotropic Luttinger model in half-space. It 
contains the usual incident and reflected waves as well as novel localized light-hole states. These surface states 
contribute to the net spin accumulation in response to a charge current. To experimentally observe the predicted 
surface states, it may be useful to consider a set-up in which a three-dimensional semiconductor (e.g., GaAs) exists 
in a close proximity to a low-dimensional system (which does not necessarily have a strong spin-orbit coupling) with 
a contact near the three-dimensional boundary. In this case, applying current to the system would result in spin 
injection into the low-dimensional system, which should be directly observable by the means of the Kerr effect. The 
experimental verification of the predicted spin-polarized surface states in hole-doped semiconductors is clearly called 
for. 
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APPENDIX A: SPIN 3/2 MATRICES AND SPINORS 



The expressions for the S matrices in a basis with the total angular momentum parallel to the z axis are: 
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For a translation invariant system the eigenvectors of the Luttinger Hamiltonian arc given by equation Q. The 
four component spinors U\(jO-k) depend on the orientation, not the magnitude, of the wave vector, so each of them 
can be uniquely described by two angles. For the purpose of this article it is enough to know their expressions in a 
particular coordinate system in which the y-component of the wave vector is zero, k v = 0. Consequently, the spinors 
will depend on one angle 0, which we choose as the angle between the z-axis and the vector — k. Explicitly we have 
for e < e c 
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For an incident angle > C there are no light-holes propagating in the z direction, but localized solutions described 
by the spinors 
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APPENDIX B: SCATTERING COEFFICIENTS FOR AN INCIDENT HEAVY-HOLE WITH HELICITY 

+3/2 



The scattering coefficients Ai and Bi in Eq. JSJ) are functions of the incident angle and can be obtained by 



imposing the boundary condition tp^ {z = 0) = 0. Explicitly we have 
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where (f)(9) characterizes the light-hole wave and is given by Eq. Q). 

In the presence of the localized modes, by imposing the boundary condition, 4 } i^ +> (z = 0) = we obtain for the 
scattering coefficients in Eq. @ 
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